I idly wonder if we can just turn this around and implement internal iteration - fold - for the indices iterators, and use that instead. It's still unlikely to generalize as well as this.

One can also see that the ndarray-rand crate has no substance beyond this line:

Array::from_shape_fn(shape, |_| dist.sample(rng))

The whole thing is a wrapper for that + a method with a default for the rng. 🙂

I don't really like the name from_shape_fn_memory_order but can't think of a better one. from_shape_fn_linear ? It's likely to cause confusion. from_shape_fn_linear_index?

I think this is a good idea, and the implementation looks good. I agree, though, that from_shape_fn_memory_order needs a better name, since the existence of a method named from_shape_fn_memory_order suggests that from_shape_fn doesn't iterate in memory order. Maybe from_shape_fn_offset to indicate the the closure's parameter is the memory offset from the first element? Alternatively, we could rename from_shape_fn to from_shape_indexed_fn and add a new from_shape_fn method that doesn't pass anything to the closure.

Now it has the name from_shape_simple_fn. It more or less sacrifices the index as parameter - now the closure is called with no parameter, so that it has an easier to understand name or explanation. We can do it either way, IMO.